#include <stdlib.h>
#include <math.h>
#include <sys/times.h>
#include "mytime.h"

#define N_DEF 200
#define TOL_DEF 0.0000001
#define MAXIT_DEF 500

/****************/
int *despl;
int *count;
/****************/
/*  No se utiliza
void CompruebaSolucion(double **a, double *x, double *b, int  n) {
  int i, j;
  double sum, dif = 0.0;

  for (i = 0; i < n; i++) {
    sum = 0.0;
    for (j = 0; j < n; j++)
      sum += a[i][j]*x[j];
    dif += (b[i] - sum)*(b[i] - sum);
  }
  dif = sqrt(dif);
  printf("||b-x||_2 = %g\n", dif);
}
*/

void print_matrix (double **a, int f, int c) {
   int i,j;
   for (i=0; i < f; i++) {
      for (j=0; j < c; j++)
         printf ("%f ", a[i][j]);
      printf ("\n");
   }
}
void print_vector (int *v, int *h,  int c) {
   int j;
   printf ("vector %d --\t", v);
   /*for (j=0; j < c; j++)
       printf ("%d ", v[j]);
   printf ("\t---- ");*/
   for (j=0; j < c; j++)
       printf ("%d ", h[j]);
   printf ("\n");
}

double jacobi_sec(double **a, int n, double *b, double *x, double tol, int maxit) {
  int i, j, it = 0;
  double sum, resi = 10.0, *new_x;

  new_x = (double *) malloc(n*sizeof(double));

  for (i = 0; i < n; i++) x[i] = b[i];
  /* while ( (it < maxit) && (resi > tol) ) { */
  while (it < maxit) {
    resi = 0.0;
    for (i = 0; i < n; i++) {
      sum = 0.0;
      for (j = 0; j < n; j++) 
        sum += a[i][j]*x[j];
      resi += fabs(sum - b[i]);
      sum += -a[i][i]*x[i];
      new_x[i] = (b[i] - sum) / a[i][i];
    }
    for (i = 0; i < n; i++) x[i] = new_x[i];
    it++;
  } 

/*  LLC_printMaster("Secuencial Iteraciones: %d, resi: %g\n", it, resi); */
  free(new_x);
  return resi;  
}

double jacobi_llc(double **a, int n, double *b, double *x, double tol, int maxit) {
  int i, j, it = 0;
  double sum, resi = 10.0, *new_x;

  new_x = (double *) malloc(n*sizeof(double));

  for (i = 0; i < n; i++) x[i] = b[i];
  /* while ( (it < maxit) && (resi > tol) ) { */
  while (it < maxit) {
    resi = 0.0;
#pragma omp target device(mpi)
#pragma omp parallel for private (sum, j) reduction (+ : resi)
#pragma llc result (&new_x[i], 1) 
    for (i = 0; i < n; i++) {
      sum = 0.0;
      for (j = 0; j < n; j++) 
        sum += a[i][j]*x[j];
      resi += fabs(sum - b[i]);
      sum += -a[i][i]*x[i];
      new_x[i] = (b[i] - sum) / a[i][i];
    }
    for (i = 0; i < n; i++) x[i] = new_x[i];
    it++;
  } 

/*  LLC_printMaster("llc Iteraciones: %d, resi: %g\n", it, resi); */
  free(new_x);

  return resi;
}


double jacobi_mpi(int miid, double **a, int nfilas, int n, double *b, double *x, double tol, int maxit) {
  int i, j, it = 0;
  double sum, resi_local, resi, *new_x;
  int fin;

  resi = 10.0;
  new_x = (double *) malloc(n*sizeof(double));

  for (i = 0; i < n; i++) x[i] = b[i];

  /* while ( (it < maxit) && (resi > tol) ) {  */
  while (it < maxit) {
    fin = despl[miid] + count[miid];
    resi_local = 0.0;
    for (i = despl[miid]; i < fin; i++) {
      sum = 0.0;
      for (j = 0; j < n; j++) 
        sum += a[i][j]*x[j];
      resi_local += fabs(sum - b[i]);
      sum += -a[i][i]*x[i];
      new_x[i] = (b[i] - sum) / a[i][i];
    }
    
    MPI_Allgatherv (&new_x[despl[miid]], count[miid], MPI_DOUBLE, x, count, despl, MPI_DOUBLE, MPI_COMM_WORLD);
    MPI_Allreduce(&resi_local, &resi, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); 
    it++;
  }


/*  LLC_printMaster("Paralelo Iteraciones %d, Residuo: %g\n", it, resi); */

  free(new_x);

  return resi;

}



int devolverCountParaHeterogeneo (int _size, int _rank, int _nproc) 
{
    int relacion;

    if (_size == 1152)
       relacion = 1;
    else 
       if (_size == 2304)
          relacion = 2;
       else 
          relacion = 4;

    
    if (_nproc == 4) 
    {
       switch (_rank) 
       {
          case 0: return (354 * relacion);
          case 1: return (297 * relacion);
          case 2: return (167 * relacion);
          case 3: return (334 * relacion);
       }
    } 
    else 
    {
       if (_nproc == 8) 
       {
          switch (_rank) 
          {
             case 0: return (180 * relacion);
             case 1: return (180 * relacion);
             case 2: return (149 * relacion);
             case 3: return (349 * relacion);
             case 4: return (79 * relacion);
             case 5: return (79 * relacion);
             case 6: return (168 * relacion);
             case 7: return (168 * relacion);
          }
       } 
       else 
       {
          switch (_rank) 
          {
             case 0: return (118 * relacion);
             case 1: return (118 * relacion);
             case 2: return (97 * relacion);
             case 3: return (97 * relacion);
             case 4: return (97 * relacion);
             case 5: return (97 * relacion);
             case 6: return (50 * relacion);
             case 7: return (50 * relacion);
             case 8: return (107 * relacion);
             case 9: return (107 * relacion);
             case 10: return (107 * relacion);
             case 11: return (107 * relacion);
          }
       }
    }
    return 0;
}

int main(int argc, char *argv[]) {
  double tiempo, tol;
  double **a, *b, *x;
  int i, j, n, nfilas, diag, maxit;
  int nprocs, miid;

  CLOCK_TYPE chrono;
  double seq_time, llc_time, mpi_time;
  double seq_resi, llc_resi, mpi_resi;

  nprocs = LLC_NUMPROCESSORS;
  miid = LLC_NAME;
 
  n = N_DEF;
  maxit = MAXIT_DEF;
  tol = TOL_DEF;

  switch (argc) {
    case 4: tol = atof(argv[3]);
    case 3: maxit = atoi(argv[2]);
    case 2: n = atoi(argv[1]);
  }

  LLC_printMaster ("\n\nEjecutando jacobi con n = %d, maxit = %d, tol = %g\n\n", n, maxit, tol);
  
  /* Por simplicidad supondremos que n es divisible por el numero de procesos: nprocs */
  nfilas = n/nprocs;

  ReservaMatriz(&a, n, n);
  ReservaVector(&b, n);
  ReservaVector(&x, n);

  if (miid == 0) {
    GeneraVector(b, n);
  }
  MPI_Bcast(b, n, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  /* matriz aleatoria */
  if (miid == 0) {
     GeneraMatriz(a, n, n);
     /* Para garantizar la convergencia del metodo de Jacobi */
     /* Convierte en diagonal dominante la matriz distribuida */
     diag = 0;
     for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) 
           a[i][diag] += a[i][j];
        diag++;
     }
  }
  /* envio de la matriz */
  for (i=0; i < n; i++) MPI_Bcast(a[i], n, MPI_DOUBLE, 0, MPI_COMM_WORLD); 
  MPI_Barrier (MPI_COMM_WORLD); 
 
  /****************************************************************************
   algoritmo secuencial */
  CLOCK_Start(chrono);
  seq_resi = jacobi_sec(a, n, b, x, tol, maxit);
  CLOCK_End(chrono, seq_time);
  MPI_Barrier (MPI_COMM_WORLD); 
 /* CompruebaSolucion (a, x, b, n); */

  /****************************************************************************
   jacobi llc */
  CLOCK_Start (chrono);
  llc_resi = jacobi_llc(a, n, b, x, tol, maxit);
  CLOCK_End(chrono, llc_time);
  MPI_Barrier (MPI_COMM_WORLD); 
/*  CompruebaSolucion (a, x, b, n);*/

  /***************************************************************************
   paralelo homogeneo */

  despl = (int *) malloc( nprocs * sizeof(int));
  count = (int *) malloc( nprocs * sizeof(int));
  despl[0] = 0;
  for (i=0; i< nprocs; i++) {
     count[i] = nfilas;
     if (i) despl[i] = despl[i-1] + count[i-1];
  }
 
  MPI_Barrier (MPI_COMM_WORLD);
  CLOCK_Start (chrono);
  mpi_resi = jacobi_mpi(miid, a, nfilas, n, b, x, tol, maxit);
  CLOCK_End(chrono, mpi_time);
  LLC_printMaster("Paralelo Procesador: %d, Size: %d, Duracion: %g\n", miid, n, tiempo);
  MPI_Barrier (MPI_COMM_WORLD); 
  /*CompruebaSolucion (a, x, b, n);*/
  
  
  /****************************************************************************
  ****************************************************************************/
  
  
  LLC_printMaster ("%d\t%g\t#llc_plot0 RAP: N=%d, IT=%d, TOL=%g (t_seq=%g/t_llc=%g) SEQ_SOL=%g, LLC_SOL=%g\n",
                            LLC_NUMPROCESSORS, seq_time/llc_time, n, maxit, tol, seq_time, llc_time, seq_resi, llc_resi);
  LLC_printMaster ("%d\t%g\t#llc_plot1 RAP: N=%d, IT=%d, TOL=%g (t_seq=%g/t_mpi=%g) SEQ_SOL=%g, MPI_SOL=%g\n",
                            LLC_NUMPROCESSORS, seq_time/mpi_time, n, maxit, tol, seq_time, mpi_time, seq_resi, mpi_resi);
  LLC_printMaster ("%d\t%g\t#llc_plot2 RAP: N=%d, IT=%d, TOL=%g (t_mpi=%g/t_llc=%g) MPI_SOL=%g, LLC_SOL=%g\n",
                            LLC_NUMPROCESSORS, mpi_time/llc_time, n, maxit, tol, mpi_time, llc_time, mpi_resi, llc_resi);


  
  /**************************************************************************
   Libeacion de memoria */
  LiberaMatriz(a);
  LiberaVector(b);
  LiberaVector(x);

  return 0;

}
